Application of the Variational Method to the Large Deformation Problem of Thin Cylindrical Shells with Different Moduli in Tension and Compression

In this study, the variational method concerning displacement components is applied to solve the large deformation problem of a thin cylindrical shell with its four sides fully fixed and under uniformly distributed loads, in which the material that constitutes the shell has a bimodular effect, in comparison to traditional materials, that is, the material will present different moduli of elasticity when it is in tension and compression. For the purpose of the use of the displacement variational method, the physical equations on the bimodular material model and the geometrical equation under large deformation are derived first. Thereafter, the total strain potential energy is expressed in terms of the displacement component, thus bringing the possibilities for the classical Ritz method. Finally, the relationship between load and central deflection is obtained, which is validated with the numerical simulation, and the jumping phenomenon of thin cylindrical shell with a bimodular effect is analyzed. The results indicate that the bimodular effect will change the stiffness of the shell, thus resulting in the corresponding change in the deformation magnitude. When the shell is relatively thin, the bimodular effect will influence the occurrence of the jumping phenomenon of the cylindrical shell.


Introduction
Thin shell structures are favored by structural engineers because of their beautiful shape and material saving. Thin shells can make full use of the strength of materials and combine the bearing and enclosure of the structure into one; in addition, under the external loads, the shell can uniformly distribute the pressure to the various parts of the shell body, making the stress state more reasonable. These advantages enable structural engineers to extensively apply thin shells in structural design [1]. At the same time, because the shell is thin, the large deformation problem often occurs; meanwhile, the stress state of thin shells is mainly compressive, prone to the stability problem. This makes the mechanical analysis of thin shells one of the topics of interest to engineers [2]. On the other hand, with the development of materials technology, thin shells composed of advanced materials are also widely used in aircraft manufacturing and ship engineering. Due to the richness and diversity of materials, the thin shell structures can serve in more demanding environments. Therefore, the research of composite thin shells has also become a hot issue in recent decades [3,4].
Still aiming for the thin shell, in this study, a new research work considering different tensile-compressive moduli (bimodular effect) of materials [5] will be conducted. Specifically, the large deformation problem of a thin cylindrical shell with a bimodular effect is proposed first and solved via the variational method. In fact, the bimodular effect of

Background
A material that has different tensile-compressive moduli is referred to as a bimodular material [9]. Most materials-for example, ceramics, concrete, rubber, graphite and some biomedical materials-exhibit different tensile-compressive strains when they are subjected to tensile-compressive stresses of the same magnitude. Generally, there are two models widely adopted in the theoretical analysis within engineering professions. The first model is proposed by Bert [10], which is established based on the criterion of positive-negative signs in the longitudinal strain of fibers. In the analysis of orthotropic materials and laminated composites [11][12][13], the Bert model has been extensively used. The second model is the Ambartsumyan model [5] established based on the criterion of positive-negative signs of principal stresses, which is applicable to isotropic materials. The Ambartsumyan model is of particular significance in structural analysis, because it is this criterion that dictates whether a certain point of structure is tensile or compressive. The present work is carried out on the Ambartsumyan model. Because the principal stress state of a point is generally obtained as a final result but not as a known condition before solving, there is no way to use the Ambartsumyan model established based on the criterion of principal stress. Analytical solutions are available in a few simple cases, although they only concern single components-for example, beams and plates [14][15][16]. In complex cases, the finite element method (FEM) on an iterative technique has to be resorted to. In each iteration, the first step is to judge the principal stress state of each element, thus obtaining a new elastic matrix used for the next iteration. This method is called the direct iterative method with variable stiffness, which was widely adopted in some earlier studies, as indicated in the reviews of Ye et al. [17] and Sun et al. [18]. Given that the traditional iteration method struggled with the convergence difficulty, Du et al. [19] successfully established a new computational framework for these kinds of constitutive models. Thereafter, based on an improved constitutive model and combined with the arc-length method, Ma et al. [20] established a finite element iterative program for bimodular rods to obtain buckling critical loads.
Among all shell structures, cylindrical shells are a special type of shells, which is characterized by zero curvature along the longitudinal direction (the bus direction of cylinder surface); this brings convenience in the analysis, design, fabrication and construction of shell structures. Thus, cylindrical shells have been regarded as an effective structural element, which is widely used in roof construction, pressure vessels and other engineering structures. There are many review works in this field, and they involve various aspects of the analysis and application of cylindrical shells, such as the earlier reviews on the buckling of moderately thick laminated cylindrical shells [21] and on the development of layered vessels using flat-ribbon-wound cylindrical shells [22], and recent reviews on the analysis of graphene nanoplatelet-reinforced cylindrical shells subjected to thermo-mechanical loads [23], on the recent progress in lightweight carbon fiber-reinforced lattice cylindrical shells [24], on the ring stiffened cylindrical shell structures [25] and on the knockdown factor of buckling load for axially compressive cylindrical shells [26], to list but a few. More recently, the natural frequencies optimization of circular cylindrical shells using axially functionally graded materials was investigated [27]. For a prestressed cylindrical shell with various structural parameters, the experimental work of the stress state was presented [28]. For sandwich shells with functionally graded coatings while operating under different external pressures, its buckling behavior was generally investigated under simply supported boundary conditions [29]. In the existing analyses of cylindrical shells, however, few works have been found on the large deformation problems of cylindrical shells, especially on the bimodular effect of materials mentioned in this study.
In general, the analytical methods of plate and shell problems fall into the following three categories: the first is the series expansion method using different function forms (for example, power function and trigonometric function); second, variational methods established based on energy principle (for example, the Ritz method and the Galerkin method); and the third is the perturbation technique. These methods have their own advantages and disadvantages which are not the focus in this paper. As solutions to the large deformation problems of plate and shell, both the variational method and the perturbation method show their unique advantages. They have the ability to provide a more consistent solution, although the approaches obtaining the same solution are quite different. In the perturbation method, the governing equation containing the unknown functions (for example, displacement and stress) should be established first. The undetermined functions are then spread in the form of ascending powers to a certain small parameter. A series of equations determining the approximate solution of all levels are thus derived by substituting the expansion into the governing equations and boundary conditions, and then by equating the same order of the perturbation parameter. Considering that the perturbation parameter either appears explicitly or is introduced artificially, Chien [30] selected the central deflection of a thin plate as a perturbation parameter to derive successfully perturbation solution in 1947. Comparing with the experimental data, Chien's solution is accurate and considered as a landmark, which has been cited in subsequent studies for a long period of time. In the variational method concerning displacement, first the displacement function containing undetermined coefficients is prescribed, which should satisfy the boundary conditions. At the same time, the functional of energy should be established, in which the total strain energy stored in the elastic body and the work performed by external loads are determined. By substituting the prescribed displacement into the functional of energy, the variation of displacement is only realized by the variation of coefficients, thus determining the unknown coefficients. In Chien's perturbation solution of the large deflection problem of a thin circular plate, the relationship between load and central deflection gives qa 4 /64D = w 0 [1 + 0.544(w 0 /t) 2 ] [30], in which q is the uniformly distributed load, a is the radius of the thin circular plate, D is the bending stiffness, t is the plate thickness and w 0 is the central deflection. For the same problem, the variational method gives the solution qa 4 /64D = w 0 [1 + 0.486(w 0 /t) 2 ] [31]. It is easy to see that they are quite close. In fact, in the previous study [32], the variational method was successfully used in the analysis of stability of cantilever vertical plates with different moduli. In addition, the recent study [33] also showed that, for the large deformation of thin shallow spherical shells, the perturbation solution agrees well with the variational solution. All these works showed that the variational method can be used to solve the large deformation problem of plate and shell.
In conclusion, the variational method has two advantages in this study. First, since the displacement variational equation itself stands for the equation of equilibrium and stress boundary conditions, it naturally avoids the equilibrium condition of shells. It is known that for cylindrical shell problems, the establishment of the equation of equilibrium is somewhat complicated. Second, the analytical characteristic of the variational method also makes it more conveniently used in the sequent parametric study and in the design of the shell. In this study, the variational method of displacement is applied to solve the large deformation problem of bimodular cylindrical shells. To this end, this paper is organized as follows. In Section 3, the variational method used and the cylindrical shell problem studied are briefly described. In Sections 4 and 5, the physical equations on the bimodular material model and the geometrical equations under the large deformation are derived, respectively. In Section 6, the total strain potential energy will be first derived and then the Ritz method is applied to solve the large deformation of bimodular cylindrical shells, thus the relationship between load and central deflection is obtained. Section 7 contains the numerical simulation and comparison with theoretical results. In Section 8, the influences of relevant parameters on this relation as well as the bimodular effect on the jumping phenomenon are discussed in detail. The concluding remarks are summarized in Section 9.

The Displacement Variational Method Based on Energy Conservation
In a three-dimensional spatial problem of elasticity whose orthogonal coordinate system is established as o-xyz, let σ x , σ y , σ z , τ xy , τ yz and τ zx be the six stress components and let ε x , ε y , ε z , γ xy , γ yz and γ zx be the six corresponding strain components, so that the strain potential energy of the whole elastic body, U, will be expressed in terms of stress and strain components as [31] Via the geometrical equation and the physical equation of the three-dimensional problem, Equation (1) may also be expressed in terms of the displacement components [31],  (2) in which u, v and w are the displacement components along the x, y and z directions, respectively. Obviously, Equation (2) opens the possibilities for the use of the displacement variational method based on energy conservation. Suppose that an arbitrary elastic body is under the action of the external forces including the body force along the x, y and z directions, X, Y and Z, and the surface force along the x, y and z directions, X, Y and Z, and now it is in equilibrium. The resulting displacement components, u, v and w, should satisfy the differential equation of equilibrium expressed in terms of displacement, displacement boundary conditions as well as stress boundary conditions expressed in terms of displacement. If the displacements take place minor changes allowed by boundary conditions, the new displacement will become where δu, δv and δw are the virtual displacement or displacement variation. During the virtual displacement, if there are no any changes in thermal energy and kinetic energy, according to the principle of conservation of energy, the increment of strain potential energy, δU, should be equal to the work done by the external forces, thus the following displacement variational equation [31] may be obtained which also may be referred to as the Lagrangian variational equation. This variational equation provides an approximate solution to problems of elasticity. Specifically, if a group of displacement components containing a series of undetermined coefficients may satisfy the displacement boundary conditions, Equation (4) may be adopted to determine these unknown coefficients, thus obtaining the final displacement. If the expression of the displacement component is taken as follows: in which A m , B m and C m are the independent coefficients whose number is 3m; u 0 , v 0 and w 0 are the specified functions whose boundary values are equal to the known displacements on the boundaries; and u m , v m and w m are specified functions which are equal to zero on the boundaries. Thus, no matter what the coefficients, A m , B m and C m , are taken, the displacement components, u, v and w, always satisfy the boundary conditions. Note that the displacement variation is realized only by the variation of A m , B m and C m , according to Equation (5), the variation of the displacement components are While the variation of the strain potential energy is Substituting Equations (6) and (7) into Equation (4), and also merging these terms with the same coefficient, Equation (4) may be changed as Since the variations, δA m , δB m and δC m , are completely arbitrary and independent each other, the coefficients of these variations in Equation (8) must be equal to zero, thus the following three relations may be obtained which may be used for the solution of the undermined coefficients, thus the displacement components may be finally obtained via Equation (5). In many studies [31,34,35], this displacement variational method is also referred to as the Ritz method.

Describtion of Problem
As shown in Figure 1, a bimodular thin cylindrical shell with its four sides fully fixed is subjected to a vertically uniformly distributed load q, in which the four corner points of the shell are expressed in points A, B, C, and D. The orthogonal curvilinear coordinate systems (θ, β, γ) is established, in which the β coordinate axis is placed along the longitudinal direction, that is, the bus direction of the cylinder surface; the θ coordinate axis (in radians) is placed along the circumferential direction, that is, the quasi-line direction of the cylinder surface; and the γ coordinate axis is placed on the direction normal to the curved plane (θ, β), whose positive direction is along the direction of the concave surface of the shell (downward). To facilitate the following analysis, the origin o is placed at the center of the shell, as shown in Figure 1. The overall size of the shell is as follows: the length of the shell along the β direction is 2b, the arc length of the shell along the θ direction is 2Ra in which a (in radians) is the half of the central angle of the cylindrical surface and R is the radius of the cylindrical surface. Note that except the load q whose unit is Pa, the units of β, γ, b and R are all length, while θ and a are in radians thus they are dimensionless quantities. Due to the geometrical characteristic of the cylindrical shell, the curvature along the θ direction is k 1 = 1/R, while the curvature along the β direction is k 2 = 0.
longitudinal direction, that is, the bus direction of the cylinder surface; the θ coordinate axis (in radians) is placed along the circumferential direction, that is, the quasi-line direction of the cylinder surface; and the γ coordinate axis is placed on the direction normal to the curved plane (θ, β), whose positive direction is along the direction of the concave surface of the shell (downward). To facilitate the following analysis, the origin o is placed at the center of the shell, as shown in Figure 1. The overall size of the shell is as follows: the length of the shell along the β direction is 2b, the arc length of the shell along the θ direction is 2Ra in which a (in radians) is the half of the central angle of the cylindrical surface and R is the radius of the cylindrical surface. Note that except the load q whose unit is Pa, the units of β, γ, b and R are all length, while θ and a are in radians thus they are dimensionless quantities. Due to the geometrical characteristic of the cylindrical shell, the curvature along the θ direction is k1 = 1/R, while the curvature along the β direction is k2 = 0. It should be noted here that the constraints for the thin cylindrical shell is considered as fully fixed on four sides. Doing so is based on the following two reasons. First, to facilitate the application of the variational method of displacement, the prescribed displacement should satisfy all displacement boundary conditions and fully fixed conditions is easily expressed in terms of the displacement components. Second, from the perspective of the application, most of the roof structures tend to be the constraints such as this, with a certain representative.
The cylindrical shell considered in this study has two distinct features, as compared to the other cylindrical shells. One is the shell will undergo the large deformation and another is the material composed the shell has a bimodular effect. From the above analysis, it is easy to see that for the application of the displacement variational method, all stress and strain components should be expressed in the form of the displacement components; thus in the next analysis, the first step is to derive the physical equation on the bimodular material model and the geometrical equation under large deformation. Due to the introduction of two features mentioned above, the physical equation and the geometrical equation will inevitably change to some extent, compared to the common equations. It should be noted here that the constraints for the thin cylindrical shell is considered as fully fixed on four sides. Doing so is based on the following two reasons. First, to facilitate the application of the variational method of displacement, the prescribed displacement should satisfy all displacement boundary conditions and fully fixed conditions is easily expressed in terms of the displacement components. Second, from the perspective of the application, most of the roof structures tend to be the constraints such as this, with a certain representative.
The cylindrical shell considered in this study has two distinct features, as compared to the other cylindrical shells. One is the shell will undergo the large deformation and another is the material composed the shell has a bimodular effect. From the above analysis, it is easy to see that for the application of the displacement variational method, all stress and strain components should be expressed in the form of the displacement components; thus in the next analysis, the first step is to derive the physical equation on the bimodular material model and the geometrical equation under large deformation. Due to the introduction of two features mentioned above, the physical equation and the geometrical equation will inevitably change to some extent, compared to the common equations.

Physical Equations on Bimodular Materials Model
The material considered in this work has a bimodular effect, that is, when the material is in tension or in compression, it will behave different Young's moduli of elasticity. As shown in Figure 1, when a thin cylindrical shell is subjected to the uniformly distributed loads acting vertically downward, there will be two different responses generated to resist the external loads, one is the bending effect and another is the in-plane compression effect, as shown in Figures 2a and 2b, respectively. These two effects act independently and not affect each other, so they can be considered separately. In other words, in each effect, if the stress state is known-for example, whether it is pulled or pressed-then the elastic modulus (tensile modulus or compressive one) can be selected accordingly. More specifically, when a certain section of the shell is under pure bending, as shown in Figure 2a, in which M is the bending moment, o' is the center of bending moment, the resulting tension-compression moduli may be determined like this. Bounded by the unknown neutral layer (which is marked in dashed line), the area above the neutral layer is in compression and the stress is compressive σ − thus the corresponding elastic modulus should be E − ; while the area below the neutral layer is in tension and the stress is tensile σ + thus the corresponding elastic modulus should be E + . Due to the tension-compression subareas under pure bending, the whole thickness of the shell, t, is divided into the tensile thickness and compressive one, which is denoted by t 1 and t 2 , respectively, as shown in Figure 2a. The determination of t 1 and t 2 may come from the bimodular thin plate under pure bending-for example, the previous study gives [36] where µ + and µ − are the tensile and compressive Poisson's ratio, respectively. In addition, Figure 2b shows the stress state of the section of the shell under in-plane loads. From the load type depicted in Figure 1, the in-plane stresses are basically compressive, thus the elastic modulus along the whole thickness may be taken as E − . modulus (tensile modulus or compressive one) can be selected accordingly. More spe ically, when a certain section of the shell is under pure bending, as shown in Figure 2a which M is the bending moment, o ' is the center of bending moment, the resulting tensio compression moduli may be determined like this. Bounded by the unknown neutral lay (which is marked in dashed line), the area above the neutral layer is in compression a the stress is compressive σ − thus the corresponding elastic modulus should be E − ; wh the area below the neutral layer is in tension and the stress is tensile σ + thus the cor sponding elastic modulus should be E + . Due to the tension-compression subareas und pure bending, the whole thickness of the shell, t, is divided into the tensile thickness a compressive one, which is denoted by t1 and t2, respectively, as shown in Figure 2a. T determination of t1 and t2 may come from the bimodular thin plate under pure bending for example, the previous study gives [36] 2 1 where μ + and μ − are the tensile and compressive Poisson's ratio, respectively. In additi Figure 2b shows the stress state of the section of the shell under in-plane loads. From load type depicted in Figure 1, the in-plane stresses are basically compressive, thus elastic modulus along the whole thickness may be taken as E − .
(a) (b) According to the basic assumptions of the common thin cylindrical shell, the tra verse shear stresses are not considered in the problem since the shell is thin and the d formation are mainly due to bending and torsion, thus ignoring the transverse sh stresses is acceptable. In addition, there are three in-plane stresses acting on the mid surface of the shell and they are the longitudinal stress, the circumferential stress and shearing stress. If let subscript 1 denote the longitudinal, let subscript 2 denote the c cumferential, let subscript 12 denote the corresponding shearing quantity between lon tudinal and circumferential directions, and the subscript t denote the in-plane, thus three in-plane stresses may be denoted by σt1, σt2 and τt12, as shown in Figure 3. Note t according to the above description, under the action of external uniformly distribut According to the basic assumptions of the common thin cylindrical shell, the transverse shear stresses are not considered in the problem since the shell is thin and the deformation are mainly due to bending and torsion, thus ignoring the transverse shear stresses is acceptable. In addition, there are three in-plane stresses acting on the middle surface of the shell and they are the longitudinal stress, the circumferential stress and the shearing stress. If let subscript 1 denote the longitudinal, let subscript 2 denote the circumferential, let subscript 12 denote the corresponding shearing quantity between longitudinal and circumferential directions, and the subscript t denote the in-plane, thus the three in-plane stresses may be denoted by σ t1 , σ t2 and τ t12 , as shown in Figure 3. Note that according to the above description, under the action of external uniformly distributed loads, the in-plane stresses acting on the whole thickness of the section tend to be compressive, thus σ t1 , σ t2 and τ t12 may be changed as σ − t1 , σ − t2 and τ − t12 and also the corresponding modulus of elasticity and Poisson's ratio may be selected as E − and µ − . The physical equation concerning the in-plane deformation may be expressed as [31] where ε t1 , ε t2 and ε t12 denote the corresponding strains.
where εt1, εt2 and εt12 denote the corresponding strains. On the other hand, the bending and torsion deformation of t be considered under large deformation. According to the abov equations concerning bending stress and torsional stress will giv of tension and compression, On the other hand, the bending and torsion deformation of the cylindrical shell must be considered under large deformation. According to the above analysis, the physical equations concerning bending stress and torsional stress will give as follows, in subareas of tension and compression, and where the subscript m denotes the out-of-plane bending or torsion; to differentiate from the subscript t (denote the in-plane), the superscript + denote the tensile and superscript − denote the compressive. σ m1 and σ m2 are the bending stresses, their resultants being the bending moments; τ m12 is the shear stress, its resultant being the torsion moment.

Geometrical Equations under Large Deformation
Generally speaking, the geometrical relations of thin cylindrical shells include two parts: one is the geometrical relation concerning bending and torsion moments, which is characterized by the second derivative of the deflection w; another involves the geometrical relation concerning in-plane deformation, not only relating to the u and v, but also to the w. However, most of the existing studies deal only with the small deformation. Under large deformation, the first geometrical relation mentioned above keeps unchanged while the second geometrical relation changes greatly, thus the geometrical equation under the large deformation is relatively complicated.

Geometrical Equation Concerning Bending and Torsion Moments
The geometrical relation concerning bending and torsion moments may be expressed in terms of the curvature due to bending (double curvature with same sign) and the curvature due to torsion (anticlastic-double curvature with opposite sign) [31] as follows, where ε m1 , ε m2 and ε m12 are the strains, w is the deflection and R is the radius of the cylindrical surface. Equation (14) may be used for the computation of the strain energies in bending and torsion.

Geometrical Equation Concerning In-Plane Deformation
Under large deformation, the strain of the middle surface generally involves the following three aspects of deformation: the first is the strain due to the displacements u and v along the directions θ and β, which is denoted by ε*; the second is the strain due to the change in the curvature radius of shell introduced by w, denoted by ε ** ; and the third one is the strain ε *** also introduced by w, but only under large deformation. That is, the final strain ε t will be the sum of the three strains, First, the middle surface strain due to the displacements u and v give [31] The strain ε 1 ** (along the θ direction) due to the change in the curvature radius of shell introduced by w may refer to Figure 4, in which the arc segment AB along the θ direction is isolated from the shell body and moves to the arc segment A 1 B 1 , after taking place a positive displacement w. The curvature radius of the shell is denoted by R, thus the initial curvature of the arc segment AB before deformation is k 1 = 1/R. The new curvature radius after deformation is R−w, and dθ is the central angle of the arc segment AB, thus the normal strain ε 1 ** along the θ direction introduced by the positive deflection w, will give Materials 2023, 16, 1686 10 of 33 the initial curvature of the arc segment AB before deformation is k1 = 1/R. The new curvature radius after deformation is R−w, and dθ is the central angle of the arc segment AB, thus the normal strain ε1 ** along the θ direction introduced by the positive deflection w, will give Figure 4. The strain ε1 ** due to the change in the curvature radius of shell introduced by w.
Due to the geometrical characteristic of the cylindrical shell, the curvature k2 along the β direction is zero, resulting in the other strain introduced by w, ε2 ** , is zero. In addi- Due to the geometrical characteristic of the cylindrical shell, the curvature k 2 along the β direction is zero, resulting in the other strain introduced by w, ε 2 ** , is zero. In addition, the shear strain, ε 12 ** , represents the change in the right angle between the θ and β directions. Now, the strain along the θ direction gives −k 1 w, which always takes place in the plane normal to the β direction while the strain along the β direction is always zero, thus the right angle between the θ and β directions cannot change, the shear strain, ε 12 ** , should be zero. Thus, The strain ε 1 *** (along the θ direction) introduced by w may refer to Figure 5, in which the arc segment AB along the θ direction whose arc length is ds, now changes to AB 1 after taking place a rotation change ∂w/∂s. Thus, the normal strain along the θ direction may be computed as Materials 2023, 16, 1686 11 of 33 Figure 5. The strain ε1 *** due to slope change in the segment along the θ direction introduced by w.
Spreading the expression in the square root by Taylor series, that is, Note that Equation (20) may be written as Similarly, ε2 *** due to the slope change in the segment along the β direction introduce by w, may be obtained as follows The shearing strain ε12 *** introduced by w may refer to Figure 6, in which the differential element ABDC changes to A1B1D1C1, after taking place the deflection w. The initial included angle between AB and AC is π/2, and after taking place the shear strain ε12 *** , the included angle between A1B1 and A1C1 is π/2−ε12 *** . Additionally, let the distance between points B and C be ds while the distance between points B1 and C1 be ds1. Thus, according to Figure 6, there is the following geometrical relation Spreading the expression in the square root by Taylor series, that is, Note that ∂w ∂s Equation (20) may be written as Similarly, ε 2 *** due to the slope change in the segment along the β direction introduce by w, may be obtained as follows The shearing strain ε 12 *** introduced by w may refer to Figure 6, in which the differential element ABDC changes to A 1 B 1 D 1 C 1 , after taking place the deflection w. The initial included angle between AB and AC is π/2, and after taking place the shear strain ε 12 *** , the included angle between A 1 B 1 and A 1 C 1 is π/2−ε 12 *** . Additionally, let the distance between points B and C be ds while the distance between points B 1 and C 1 be ds 1 . Thus, according to Figure 6, there is the following geometrical relation In the triangle A 1 B 1 C 1 , the length of ds 1 may be computed as follows, via cosine theorem, where and similarly, Spreading the expression in the square root by Taylor series, that is, Note that Equation (20) may be written as Similarly, ε2 *** due to the slope change in the segment along the β direction introduce by w, may be obtained as follows The shearing strain ε12 *** introduced by w may refer to Figure 6, in which the differential element ABDC changes to A1B1D1C1, after taking place the deflection w. The initial included angle between AB and AC is π/2, and after taking place the shear strain ε12 *** , the included angle between A1B1 and A1C1 is π/2−ε12 *** . Additionally, let the distance between points B and C be ds while the distance between points B1 and C1 be ds1. Thus, according to Figure 6, there is the following geometrical relation Figure 6. The shearing strain ε12 *** introduced by w. Figure 6. The shearing strain ε 12 *** introduced by w.
Combining Equations (24) and (25) yields Due to 1 R ∂w ∂θ the following strain will be obtained Thus, collecting the strain only under large deformation introduced by w will yield By adding Equations (16), (18) and (31), the geometrical equation concerning middle surface deformation is finally obtained as follows

Total Strain Potencial Energy
The total strain potential energy, U, consists of the strain energy due to the deformation of middle surface, U t , and the energy due to the deformation of bending and torsion, U m , that is [31], where the subscript t denotes the in-plane and the subscript m denotes the bending or torsion, corresponding to the notational conventions in the geometrical and physical equations above. First, the strain energy due to the deformation of the middle surface, U t , is computed as follows [31] Substituting Equation (11) into Equation (34) and also noting that the lower limit to upper limit of integration along the γ direction gives -t 2 to t 1 , Equation (34) may be written as Additionally, substituting Equation (32) into Equation (35) will yield Thus, the U t is expressed in terms of u, v and w.
The energy due to the deformation of bending and torsion, U m , may be derived, via the subareas in tension and compression indicated above, that is, Substituting Equations (12)-(14) into the above equation, also noting that integration limits in the tensile term is from 0 to t 1 while the limits in the compressive term is from −t 2 to 0, will yield Integrating with respect to γ will yield Due to dS = Rdθdβ, and according to the Green formula, Equation (40) is changed as Because the four sides of the cylindrical shell are fully fixed, thus yielding Let which is the bending stiffness of a bimodular cylindrical shell in bending; also let the Laplace operator be [31] the strain energy U m may be finally expressed as Combining Equations (36) and (46), the total strain potential energy which is expressed in terms of the displacement components is finally obtained.

The Ritz Method
For the large deformation problem of the cylindrical shell, the expression of the displacement components is taken as follows: in which A m , B m and C m are the independent coefficients; and u m , v m and w m are specified functions which are equal to zero on the boundaries. Thus, the displacement components, u, v and w, always satisfy boundary conditions of displacement. As indicated above, the four sides of the cylindrical shell are fully fixed, thus the boundary conditions of displacement give where a (in radians) is the angle along the θ direction and b is the length along the β direction, please see Figure 1. Further, the displacement u and v also satisfy the symmetry conditions, that is, this can be easily observed from Figure 1. According to the analysis, the specific expressions of the displacement components can be supposed as follows and Equations (50)-(52) can satisfy all the boundary conditions of displacement, that is, Equations (48) and (49). Specifically, the term θ 2 −a 2 can satisfy u, v and w are zero when θ = ±a; similarly, the term β 2 −b 2 can satisfy u, v and w are zero when β = ±b; the single θ and β in u and v are used for satisfying the conditions that u and v are zero when θ and β are zero, respectively. The last term in Equations (50)-(52), which is expressed in terms of the series of θ 2 and β 2 , is mainly used for the asymptotic in mathematics. In addition, u, v and w are the even functions or odd functions with respect to θ and β, this feature is also easily observed from Equations (50)-(52).
According to the conclusion from [31,34,35], taking only the first few terms, the variational method of displacement can give relatively satisfactory results. In the next computation, for convenience, the first three undetermined coefficients, A 0 , A 1 and A 2 , in Equation (50) may be taken; this is the same practice for the three coefficients, B 0 , B 1 and B 2 , in Equation (51) while the first coefficient C 0 in Equation (52) is taken; and then substituting these displacement formulas into Equations (36) and (46) to obtain the total strain potential energy, which is expressed in terms of A 0 , A 1 , A 2 , B 0 , B 1 , B 2 and C 0 .
According to the variational method of displacement described in Section 3.1 and also considering Equation (9), the following three variational equations will be obtained and where Z = q. Via Equations (53) and (54), the equations used for the solution of the coefficients A 0 , A 1 , A 2 , B 0 , B 1 , and B 2 may be obtained as follows, By using the mathematical software Maple 2021, the coefficients A 0 , A 1 , A 2 , B 0 , B 1 and B 2 may be computed, all of which contain C 0 and are listed in Appendix A.
At the same time, when C 0 is considered only, Equation (52) is changed as Substituting it into Equation (55) and integrating will yield Because A 0 , A 1 , A 2 , B 0 , B 1 and B 2 are now known, substituting them into Equation (58), an equation which contains only C 0 may be obtained. On the other hand, according to Equation (52), and also let the central deflection, that is, the maximum deflection, be w 0 when θ and β are zero, Equation (52) will give This means that the central deflection w 0 is indeed the coefficient C 0 , only with the difference of the factor a 4 b 4 . Substituting Equation (59) into Equation (58), the relation of the load vs. central deflection may be finally obtained as follows, in which H 1 (λ, µ − ) to H 5 (λ, µ − ) are the formulas containing λ and µ − , and they are listed in Appendix B, and λ = Ra/b.

Numerical Simulation and Comparison with Theoretical Solution
The numerical simulation is conducted by using the software ABAQUS6.14.4, in which the subroutine UMAT is adopted since, in the existing commercial software, there is no bimodular material model. The whole iteration process will proceed in the following steps: (i) First, the material is assumed to have a single modulus, and the stress and strain of each element are calculated; (ii) The principal stress and its direction for each element are thus obtained; (iii) According to principal stress obtained, the due constitutive relationship between each element is constructed, and the stiffness matrix of each element is collected to form the total stiffness matrix; (iv) According to the new constitutive relationship, once again, the stress and strain of each element is calculated, and repeat this process; (v) If reaching the stopping condition, output the final result; otherwise go to (ii).
When constructing the computational model of cylindrical shell, the first step is to form the three-dimensional diagram according to the real shape and size. In this study, three groups of model of thin cylindrical shells are used, in which the circumferential length of the shell, 2Ra, is 20 m, the longitudinal length, 2b, is 40 m and the thickness of the shell, t, is 0.2 m (please refer to Figure 1), with different Ra 2 /t values ranging from 1/2, 1 to 3/2. Figure 7 shows the cylindrical shell model. Note that since the cylindrical shell is relatively shallow, the side elevation and plan of the model, but the three-dimensional diagram, are given in Figure 7.
length of the shell, 2Ra, is 20 m, the longitudinal length, 2b, is 40 m and the thickness of the shell, t, is 0.2 m (please refer to Figure 1), with different Ra 2 /t values ranging from 1/2, 1 to 3/2. Figure 7 shows the cylindrical shell model. Note that since the cylindrical shell is relatively shallow, the side elevation and plan of the model, but the three-dimensional diagram, are given in Figure 7. In the numerical simulation, the average modulus, E, is taken as 20 Gpa; and the average Poisson's ratio, μ, is taken as 0.3. Five different ratios of E + /E − are selected, and they are 1/2, 2/3, 1, 3/2 and 2 (please refer to the next section). The load intensity ranges from 15 KPa to 30 Kpa, with an interval of 5 Kpa. To compare with the theoretical solution, the boundary conditions in the numerical simulation are also taken as fully fixed. A threedimensional solid element with 8 nodes, C3D8, is adopted. Some representative displacement nephograms are shown in Figures 8-10, in which Figure 8 shows the displacement nephograms when Ra 2 /t = 1/2 and q = 20 Kpa; Figure 9 shows the nephograms when Ra 2 /t = 1 and q = 25 Kpa and Figure 10 shows the nephograms when Ra 2 /t = 3/2 and q = 30 Kpa. From Figures 8-10, it is readily apparent that the maximum deflection takes place at the central part of the cylindrical shell; the deflection is gradually reduced as the observation point approaches the shell edge.  In the numerical simulation, the average modulus, E, is taken as 20 Gpa; and the average Poisson's ratio, µ, is taken as 0.3. Five different ratios of E + /E − are selected, and they are 1/2, 2/3, 1, 3/2 and 2 (please refer to the next section). The load intensity ranges from 15 KPa to 30 Kpa, with an interval of 5 Kpa. To compare with the theoretical solution, the boundary conditions in the numerical simulation are also taken as fully fixed. A three-dimensional solid element with 8 nodes, C3D8, is adopted. Some representative displacement nephograms are shown in Figures 8-10, in which Figure 8 shows the displacement nephograms when Ra 2 /t = 1/2 and q = 20 Kpa; Figure 9 shows the nephograms when Ra 2 /t = 1 and q = 25 Kpa and Figure 10 shows the nephograms when Ra 2 /t = 3/2 and q = 30 Kpa. From Figures 8-10, it is readily apparent that the maximum deflection takes place at the central part of the cylindrical shell; the deflection is gradually reduced as the observation point approaches the shell edge.
length of the shell, 2Ra, is 20 m, the longitudinal length, 2b, is 40 m and the thickness of the shell, t, is 0.2 m (please refer to Figure 1), with different Ra 2 /t values ranging from 1/2, 1 to 3/2. Figure 7 shows the cylindrical shell model. Note that since the cylindrical shell is relatively shallow, the side elevation and plan of the model, but the three-dimensional diagram, are given in Figure 7. In the numerical simulation, the average modulus, E, is taken as 20 Gpa; and the average Poisson's ratio, μ, is taken as 0.3. Five different ratios of E + /E − are selected, and they are 1/2, 2/3, 1, 3/2 and 2 (please refer to the next section). The load intensity ranges from 15 KPa to 30 Kpa, with an interval of 5 Kpa. To compare with the theoretical solution, the boundary conditions in the numerical simulation are also taken as fully fixed. A threedimensional solid element with 8 nodes, C3D8, is adopted. Some representative displacement nephograms are shown in Figures 8-10, in which Figure 8 shows the displacement nephograms when Ra 2 /t = 1/2 and q = 20 Kpa; Figure 9 shows the nephograms when Ra 2 /t = 1 and q = 25 Kpa and Figure 10 shows the nephograms when Ra 2 /t = 3/2 and q = 30 Kpa. From Figures 8-10, it is readily apparent that the maximum deflection takes place at the central part of the cylindrical shell; the deflection is gradually reduced as the observation point approaches the shell edge.  length of the shell, 2Ra, is 20 m, the longitudinal length, 2b, is 40 m and the thickness of the shell, t, is 0.2 m (please refer to Figure 1), with different Ra 2 /t values ranging from 1/2, 1 to 3/2. Figure 7 shows the cylindrical shell model. Note that since the cylindrical shell is relatively shallow, the side elevation and plan of the model, but the three-dimensional diagram, are given in Figure 7. In the numerical simulation, the average modulus, E, is taken as 20 Gpa; and the average Poisson's ratio, μ, is taken as 0.3. Five different ratios of E + /E − are selected, and they are 1/2, 2/3, 1, 3/2 and 2 (please refer to the next section). The load intensity ranges from 15 KPa to 30 Kpa, with an interval of 5 Kpa. To compare with the theoretical solution, the boundary conditions in the numerical simulation are also taken as fully fixed. A threedimensional solid element with 8 nodes, C3D8, is adopted. Some representative displacement nephograms are shown in Figures 8-10, in which Figure 8 shows the displacement nephograms when Ra 2 /t = 1/2 and q = 20 Kpa; Figure 9 shows the nephograms when Ra 2 /t = 1 and q = 25 Kpa and Figure 10 shows the nephograms when Ra 2 /t = 3/2 and q = 30 Kpa. From Figures 8-10, it is readily apparent that the maximum deflection takes place at the central part of the cylindrical shell; the deflection is gradually reduced as the observation point approaches the shell edge.     Tables 1-3 list the central deflection under different bimodular cases (E + /E − = 1/2, 2/3, 1, 3/2 and 2), different loads (q = 15 Kpa, 20 Kpa, 25 Kpa and 30 Kpa) and different shapes (Ra 2 /t = 1/2, 1 and 3/2). The results of the variational method are from Equation (60) and the FEM results are from the numerical simulation. During the mesh division, the shell is divided as 5 or 4 layers along the direction of thickness, with the difference of grid sizes. Specifically, for the case Ra 2 /t = 1/2, the grid size may be 0.2 m or 0.35 m, thus resulting in 100,000 elements or 32,490 elements, respectively; for the case Ra 2 /t = 1, the grid size may be 0.25 m or 0.35 m, thus resulting in 64,000 elements or 25,993 elements, respectively; for the case Ra 2 /t = 3/2, the grid size may be 0.24 m or 0.33 m, thus resulting in 69,305 elements or 29,524 elements, respectively.
By comparing the values of central deflection in Tables 1-3, it is easy to see that the values from two different methods are basically consistent but there still exist the differences between them. The variational solution is based on the simplified mechanical model on tension-compression subareas while the numerical simulation is based on the original mechanical model concerning positive-negative signs of principal stresses. Generally, the differences are acceptable, which validates the variational method from the side.

Load vs. Central Deflection
For the convenience of the next analyses, the following quantities are introduced where η is the ratio of the difference of tensile-compressive moduli to the sum of tensilecompressive moduli, and it is obvious that η is a small dimensionless quantity which corresponds to the positive (E + > E − ) or the negative (E + < E − ); µ and E is the average Poisson's ratio and average modulus, respectively; P and W 0 is the dimensionless load and central deflection, respectively; K is the dimensionless bending stiffness, and the relation among K, D* and A is also shown in Equation (61), and A is dimensionless. Thus, Equation (60) now becomes, noting that H 1 (λ, µ − ) to H 5 (λ, µ − ) are dimensionless in themselves, In the bimodular problem, it is generally assumed that the tensile-compressive moduli and tensile-compressive Poisson's ratios satisfy the familiar relation, E + /E − = µ + /µ − [5]. Therefore, by combining Equation (61), the tensile-compressive moduli and Poisson's ratios may be computed as In addition, the dimensionless tensile-compressive thickness, t 1 /t and t 2 /t, are determined via Equation (10) from the previous study [34]. If the η value is given in advance, E + and E − may be expressed in the form of E via Equation (63); at the same time, if µ is prescribed as 0.3, which is a representative value for most engineering materials (for example, metal), µ + and µ − may be determined; thus, t 1 /t and t 2 /t, as well as A and K may also be determined via Equation (61). These given values and the corresponding computed values are listed in Table 4.  According to Equation (62) and Table 4, the curves of P-W 0 under different λ = Ra/b values (1/2, 1 and 2), different E + /E − values (1/2, 2/3, 1, 3/2 and 2) and different Ra 2 /t values (ranging from 3 to 11/2, with the interval 1/2) may be plotted, as shown in Figures 11-13, in which Figure 11 shows the case of λ = 1/2, and Figure 12 shows the case of λ = 1, and Figure 13 shows the case of λ = 2. Figures 11-13 all show that, under the same magnitude of P, among the five cases of different moduli, the produced deflection values are in turn, from the least to the most, E + /E − = 2, 3/2, 1, 2/3 and 1/2, indicating that when E + > E − , the central deformation of the cylindrical shell will decrease and when E + < E − , the central deformation of the cylindrical shell will increase. It is obvious that the bimodular effect will change the structural stiffness of the shell, thus resulting in the corresponding change in deformation magnitudes. According to Equation (62) and Table 4, the curves of P-W0 under different λ = Ra/b values (1/2, 1 and 2), different E + /E − values (1/2, 2/3, 1, 3/2 and 2) and different Ra 2 /t values (ranging from 3 to 11/2, with the interval 1/2) may be plotted, as shown in Figures 11-13, in which Figure 11 shows the case of λ = 1/2, and Figure 12 shows the case of λ = 1, and Figure 13 shows the case of λ = 2.  For shell structures, the jumping phenomenon is prone to occur in real applications, especially for thin and shallow shells, and the bimodular cylindrical shell presented in this study is no exception. From the perspective of physical phenomena, jumping is a sudden change from one equilibrium state to another, which involves the nonlinear problem and must be solved by nonlinear solving methods. The generation of jumping is influenced by multiple factors, such as initial deflection and uneven stress distribution of shells, but this is not the focus of this study. From Figure 11, the jumping phenomenon is easily observed firstly, which is characterized by the occurrence of inflection point of the curve. However, From Figure 12, the inflection point can be found only in some cases. Further, for Figure 13, no inflection point can be found. This phenomenon indicates that the occurrence of inflection point will depend on several factors-for example, λ = Ra/b values, different E + /E − values and different Ra 2 /t values. More in-depth analysis is still needed, especially for the bimodular effect on the jumping phenomenon.   Figures 11-13 all show that, under the same magnitude of P, among the five cases of different moduli, the produced deflection values are in turn, from the least to the most, E + /E − = 2, 3/2, 1, 2/3 and 1/2, indicating that when E + > E − , the central deformation of the cylindrical shell will decrease and when E + < E − , the central deformation of the cylindrical shell will increase. It is obvious that the bimodular effect will change the structural stiffness of the shell, thus resulting in the corresponding change in deformation magnitudes.
For shell structures, the jumping phenomenon is prone to occur in real applications, especially for thin and shallow shells, and the bimodular cylindrical shell presented in this study is no exception. From the perspective of physical phenomena, jumping is a sudden change from one equilibrium state to another, which involves the nonlinear problem and must be solved by nonlinear solving methods. The generation of jumping is influenced by multiple factors, such as initial deflection and uneven stress distribution of

Jumping Phenomenon
The attention will be focused on the influence of a bimodular effect on the jumping phenomenon. To this end, the case λ = 1/2 is selected as the studied object since in this case, the jumping phenomenon is easily observed. The curves of load vs. central deflection under different Ra 2 /t values are plotted, as shown in Figures 14-16, in which Figure 14 is for the case of E + /E − = 1/2, Figure 15 is for E + /E − = 1, and Figure 16 is for E + /E − = 2.    Figure 14 shows that the jumping phenomenon may occur when Ra 2 /t is approximately 4, while from Figure 15, the Ra 2 /t value is approximately 5 and from Figure 16, the Ra 2 /t value is approximately 6, indicating that under the same λ value, the bimodular effect will change the occurrence values of Ra 2 /t of jumping phenomenon. The smaller Ra 2 /t    Figure 14 shows that the jumping phenomenon may occur when Ra 2 /t is approximately 4, while from Figure 15, the Ra 2 /t value is approximately 5 and from Figure 16, the Ra 2 /t value is approximately 6, indicating that under the same λ value, the bimodular effect will change the occurrence values of Ra 2 /t of jumping phenomenon. The smaller Ra 2 /t    Figure 14 shows that the jumping phenomenon may occur when Ra 2 /t is approximately 4, while from Figure 15, the Ra 2 /t value is approximately 5 and from Figure 16, the Ra 2 /t value is approximately 6, indicating that under the same λ value, the bimodular effect will change the occurrence values of Ra 2 /t of jumping phenomenon. The smaller Ra 2 /t  Figure 14 shows that the jumping phenomenon may occur when Ra 2 /t is approximately 4, while from Figure 15, the Ra 2 /t value is approximately 5 and from Figure 16, the Ra 2 /t value is approximately 6, indicating that under the same λ value, the bimodular Via Equation (66), the critical value of Ra 2 /t under different moduli may be calculated. When E + /E − = 1/2, the critical value of Ra 2 /t is 3.80; when E + /E − = 1, the critical value of Ra 2 /t gives 4.89 and when E + /E − = 2, the critical value of Ra 2 /t changes to be 6.08. These three theoretical values agree with the results observed from Figures 14-16.
Substituting the two real roots of Equation (64) into Equation (62), the relation of critical load Pcr and the parameter Ra 2 /t when λ = 1/2 may be obtained, which is shown in Figure 17. It is easy to see that from Figure 17, the three groups of curves under different moduli have their own sharp corners, in which the case of E + /E − = 1/2 corresponds to Ra 2 /t = 3.80, the case of E + /E − = 1 to Ra 2 /t = 4.89 and the case of E + /E − = 2 to Ra 2 /t = 6.08.

Concluding Remarks
In this study, the displacement variational method is applied to solve the large deformation problem of bimodular cylindrical shells. For the convenience of the use of the variational method, the physical equations on the bimodular material model and the geometrical equation under large deformation are derived first; the total strain potential energy is expressed in terms of the displacement component, thus bringing the possibilities for the realization of the Ritz method. Finally, the analytical relationship between load

Concluding Remarks
In this study, the displacement variational method is applied to solve the large deformation problem of bimodular cylindrical shells. For the convenience of the use of the variational method, the physical equations on the bimodular material model and the geometrical equation under large deformation are derived first; the total strain potential energy is expressed in terms of the displacement component, thus bringing the possibilities for the realization of the Ritz method. Finally, the analytical relationship between load and central deflection is obtained and validated with the numerical simulation. The jumping phenomenon of thin cylindrical shells with a bimodular effect is also discussed. The following three conclusions can be drawn.
(i) Under large deformation and using bimodular materials, the establishment of physical equations and the geometrical equation lays the foundation for the calculation of strain energy and then the application of the variational method. (ii) The bimodular effect will change the structural stiffness, thus resulting in the corresponding change in the deformation magnitude. Specifically, under the action of the same load, the central deflection of the cylindrical shell will decrease when E + > E − while the central deflection of the cylindrical shell will increase when E + < E − , in comparison with the central deflection when E + = E − . (iii) When the shell is relatively thin, the bimodular effect will influence the occurrence of the jumping phenomenon of the cylindrical shell. Specifically, if λ = 1/2, that is, the arc length along the θ direction is half of line length along the β direction, the shell is a long shell, when E + /E − = 1/2, the critical value of Ra 2 /t is 3.80; when E + /E − = 1, the critical value of Ra 2 /t gives 4.89 and when E + /E − = 2, the critical value of Ra 2 /t changes to be 6.08.
In addition, the occurrence of the jumping phenomenon depends on many factors, including λ = Ra/b values (is the cylindrical shell relatively long or short?), different Ra 2 /t values (is the cylindrical shell relatively thick or thin?), different load types (uniformly distributed loads or concentrated force, or both) and all kinds of boundary conditions (fully fixed, simply supported or the mix). Limited to the use of the variational method, the bimodular cylindrical shell in this study is fully fixed at its four sides, and the load considered is also uniformly distributed load. In practical problems, the load and boundary conditions may be diverse. More investigations are still needed, especially for the relatively thin cylindrical shells when the above factors are superimposed together.
If shell structures are composed of the materials with a relatively obvious bimodular effect; or even if the material has no obvious bimodular effect, the analytical results obtained in this study may be used for the refined analysis and optimized design. Specifically, via the analytical relationship between load and central deflection, it is readily known for the designers what will be the maximum deflection under a certain design load, which is often required in the design codes and specifications.
The limitations of this study exist in the following two aspects: the cylindrical shell is a particular case of a shell with curvature in both directions, so this general case should be studied in future studies; this study does not incorporate the effects of a snap-through phenomenon due to instability, so this effect should also be incorporated in future studies.